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Preface 


The  subject  of  low-thrust,  ion  rockets  has  interested 
ae  since  I  first  learned  of  such  vehicles  in  an  under¬ 
graduate  astrodynamics  course.  Air  Force  duties  prevented 
ae  froa  keeping  abreast  of  new  developments  over  the  past 
few  years,  however.  I  welcomed  the  opportunity  of  working 
on  this  subject  for  my  thesis. 

I  would  like  to  acknowledge  the  help  of  my  advisor, 
Dr.  William  Wiesel.  But  I  especially  want  to  thank  my 
wife  and  typist,  Jan,  for  her  understanding  and  help  over 
the  past  18  months.  I  also  hope  my  two  daughters,  Laura 
and  Leslie,  will  forgive  me  for  the  lack  of  attention  they 
have  had  to  experience. 


Accession  For 

NTIS  GRA4I 
DTIC  TAB 
Unannounced  □ 

Justification - 


By - 

Distribution/ 


Availability  Codes_ 
Avail  and/or 
Special 


Bob  Cass 


ii 


Table  of  Contents 


Preface  .  . 

List  of  Figures  . . . 

List  of  Symbols  . . . 

Abstract  .  ........ 

I.  Introduction  •  •  •  •  .  . 

II.  Fast  Timescale  Problem  . 

Derivation  . 

Examples  . 

III.  Slow  Timescale  Problem  •  •  .  • 

Problem  Statement  . 

Derivation  .  • 

Implementation  . 

IV.  Results . . . 

Case  I  . 

Case  II  . 

Case  III  ..••••••• 

Multiple  Local  Minima  •  •  • 
General  Observations  •  .  • 

V.  Conclusions  and  Recommendations  • 
Bibliography 


Figure  Page 

1 •  Shadow  Boundaries  .  7 

2.  Acceleration  Components  .  10 

3.  Sun-Orbit  Geometrical  Relationship  .  13 

4.  Shadow  Geometry . 16 

5.  Shadow  Effects  on  . 19 

6*  Shadow  Effects  on  Aa  and  Ai . 20 

7.  OL  and  y  for  m  =  0 . 23 

8.  a  and  y  for  m  =  7T/6 . 24 

9.  Oi  and  y  for  m  =  7T/3 . 25 

10.  X  a  and  X±  for  Initial  Conditions . 40 

11.  Xa  and  Xi  for  Final  Conditions . 40 

12.  Case  I  •••••••••••.••• . 43 

13.  Case  II  . 44 

14.  Case  III . 45 

15.  Initial  Conditions  Case  II . 46 

16.  Initial  Conditions  Case  III . 47 

17.  T  versus  Xa(0)  49 


List  of  Symbols 


a  semi-major  axis 

af  final  a 

a^  initial  a 

B  second  position  angle  of  anti-sun  point 
D  function  of  a,  f,  s 

DU  earth  distance  unit 

e  orbital  eccentricity 

f  true  anomaly 

F  integrand  of  performance  index 

F0  constant  in  computation  of  (J0 

g  true  anomaly  at  shadow  entry 

G  universal  gravitation  constant 

G0  constant  in  computation  of  (J0 

H  Hamiltonian 

i  orbital  inclination 

if  final  i 

J  performance  index 

L  first  position  angle  of  anti-sun  point 
m  one-half  actual  shadow  angle 

mp  specific  mass  flow  rate 

M#  mass  of  earth 

n  mean  motion 

Hq  angular  semi-diameter  of  sun 

s  angle  between  ascending  node  and  shadow  exit 
t  time 


v 


T 


total  acceleration  from  thrust 
initial  T 


earth  time  unit 

control  variable 

radial  acceleration  component 

tangential  acceleration  component 

normal  acceleration  component 

da/  da 

dm/  di 

orbital  period  of  satellite 

thrust  angle  in  orbit  plane 

right  ascension  of  sun 

thrust  angle  out  of  orbit  plane 

declination  of  sun 

Lagrange  multiplier 

gravitational  parameter 

parallax  of  satellite 

parallex  of  sun 

one-half  maximum  shadow  angle 

independent  variable  in  slow  timescale  problem 

argument  of  perigee 

longitude  of  the  ascending  node 


AFIT/GA/AA/83D-1 


Abstract 

This  paper  examines  the  use  of  discontinuous  low 
thrust  for  orbital  transfers  between  two  non-coplanar, 
circular  orbits.  The  vehicle  is  assumed  to  be  a  solar- 
powered,  ion  rocket  that  cannot  operate  when  it  is  within 
the  earth's  shadow.  Two  timescales  are  used  to  derive 
a  minimum  fuel  trajectory.  The  fast  timescale  solution 
maximizes  a  change  in  inclination. when  given  a  change  in 
semi-major  axis  for  a  single  orbit.  The  slow  timescale 
solution  combines  fast  timescale  results  to  obtain  the 
minimum  fuel  trajectory.  Results  are  presented  for 
three  specific  transfers  requiring  varying  amounts  of 
shadow  penetration.  It  is  shown  that  the  fuel  penalty 
caused  by  discontinuous  thrust  is  very  small.  However, 
there  can  be  a  moderate  increase  in  total  trip  time  if 
the  time  within  shadow  is  large. 


DISCONTINUOUS  LOW  THRUST  ORBIT  TRANSFER 

I  Introduction 

Electric  propulsion  research  continues  to  seek  new 

methods  accomplishing  future  missions  in  space. ^  Elec- 
* 

trically  propelled  vehicles  offer  two  important  advantages 
for  anticipated  missions  dealing  with  large  space  struc¬ 
tures.  One  advantage  is  an  increase  in  payload  ratio 
which  will  allow  more  massive  structures  to  be  propelled. 

The  other  advantage  comes  from  the  extremely  low  accelera¬ 
tion  available  from  electric  thrusters.  The  large  space 
structures  being  proposed  will  be  flimsy  and  unable  to 
withstand  large  accelerations.  Electric  propulsion  offers 
a  way  of  moving  these  large  structures  around  in  space. 

Although  there  are  many  versions  of  electric  thrusters, 
they  all  have  common  characteristics.  The  most  notable  are 
high  specific  impulse,  low  mass  flow  rate  and  low  thrust. 
These  devices  provide  thrust  by  electrically  accelerating 
charged  ions  and  then  exhausting  them  into  space. 

Several  authors  (Alfano^,  Moeckel^,  Ehrike\  and 
Stuhlinger^)  have  shown  that  the  optimal  trajectory  for 
orbital  transfers  using  continuous  low  thrust  is  an  out¬ 
ward  spiral*  Alfano  has  derived  an  optimal  control  law 
to  perform  transfers  that  include  changes  in  both  inclination 
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and  semi-major  axis.  All  these  solutions  have  been  based 
on  the  use  of  continuous  low  thrust. 

But  continuous  thrust  may  not  be  available.  The  electric 
power  used  to  drive  these  thrusters  will  probably  be  produced 
by  solar  cells  or  a  nuclear  power  plant.  Since  nuclear  power 
sources  weigh  more  per  kilowatt  than  solar  cells  and  since 
nuclear  power  is  becoming  less  attractive  in  general,  it  is 
reasonable  to  assume  that  many  of  these  vehicles  will  be 
solar- powered.^  A  solar-powered  rocket  has  one  serious  dis¬ 
advantage,  however.  It  will  not  work  when  it  is  in  the  earth's 
shadow.  Therefore,  such  a  craft  could  not  provide  continuous 
thrust. 

If  tangential  thrust  were  applied  only  when  the  vehicle 
is  in  sunlight,  a  circular  orbit  would  not  remain  circular 
very  long.  Since  the  previous  control  schemes  all  assumed 
circular  orbits  would  remain  circular,  they  are  not  valid 
for  discontinuous  low  thrust.  This  paper  addresses  the 
problem. of  using  discontinuous  low  thrust  for  non-coplanar 
transfers  between  circular  orbits. 

The  derivation  is  divided  into  two  problems  of  differing 
time  scales.  The  fast  timescale  problem  optimizes  the 
changes  in  orbital  elements  over  one  orbit  with  vehicle  mass 
held  constant.  The  slow  timescale  problem  uses  the  results 
of  the  fast  timescale  optimization,  while  updating  the  mass 
and  acceleration  on  each  orbit. 
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II  Fast  Timescale  Problem 


The  fast  timescale  problem  addresses  the  changes  in 
orbital  elements  that  occur  during  one  revolution  of  the 
central  body.  The  solution  should  produce  a  thrust  profile 
which  will  maximize  a  change  in  either  inclination  or  semi- 
major  axis  when  a  particular  change  in  the  other  is  speci¬ 
fied.  Thrust  application  is  only  possible  when  the  solar- 
powered  spacecraft  is  in  sunlight. 


Derivation 

The  following  derivation  assumes  two-body  motion  with 
the  earth  as  the  central  body.  Electric  engines  produce 
low  thrust  and  hence  cause  only  small  changes  in  the  orbital 
elements  during  one  orbit.  Consequently,  general  perturba¬ 
tion  theory  is  used  in  that  all  the  orbit  elements  are  con¬ 
sidered  constant  during  each  orbit.  Also,  fuel  consumption 
is  low  enough  that  mass  is  considered  constant  during  one 
orbit. 

The  equations  of  motion  are  given  by  Lagrange's  planetary 
equations  in  their  acceleration  component  form: ^ 


dQ  W  (1-e2)*  sinCf+GJ) 
dt  n  a  (1  ♦  e  cos  f)  sin  i 

di  W  (1-e2)*cos(f+0)) 
dt  n  a  (1  +  e  cos  f) 


\"  V  *.  \  »*.  * 

•  •  •  •  •-  •  -  *  »  •  •  • 


U  (l-e^)*  cos  f  V  (1 -e^)T(2  +  e  cos  f)  sin  f 


dt  nae  n  a  e  (1  e  cos  f) 

W  (1-e2)*  sin(f+(J)  cot  i 
n  a  (1  +  e  cos  f) 


(3) 


de  U  (1-e2)*  sin  f 
dt  n  a 


V  (1-e2)* 
nae 


1  -  e 


+  e  cos  f  - 


1  +  e  cos  f 


(4) 


da  2  U  e  sin  f  2  V  (1  +  e  cos  f) 

dt  n  (1-e2)^  n  (1-e2)^ 


(3) 


where  a  is  the  semi-major  axis;  e  is  the  eccentricity;  i  is 
inclination; (J  is  the  argument  of  perigee;  Q  is  the  longitude 
of  the  ascending  node;  f  is  the  true  anomaly;  U,  V  and  W  are 
the  radial,  tangential  and  normal  acceleration  components, 
respectively;  and  n  is  the  mean  motion. 

n  =  (yu/a3)*  (6) 

where 

M  =  G  Ma  (7) 


and  G  is  the  universal  gravitation  constant  and  is  the 


mass  of  the  earth* 

For  transfers  between  circular  orbits  with  continuous 
low  thrust,  many  authors  have  shown  that  the  optimum  thrust 
profile  requires  negligible  radial  thrust.  It  has  also  been 
shown  that  the  use  of  purely  tangential  thrust  causes  only 
negligible  changes  in  eccentricity  for  transfers  to  geosyn¬ 
chronous  altitude  when  the  thrust  acceleration  is  less  than 
10“^g's,^  But  when  thrust  is  applied  over  only  part  of  the 
orbit,  eventual  changes  in  eccentricity  can  be  expected.  It 
is  also  possible  that  the  optimum  thrust  profile  would  re¬ 
quire  that  eccentricity  be  allowed  to  vary  from  zero.  But 
the  planetary  equations  can  be  greatly  simplified  if  the 
dependence  on  eccentricity  can  be  eliminated. 

To  eliminate  eccentricity  from  the  planetary  equations, 
an  additional  constraint  is  added  to  the  problem.  It  will 
be  required  that  the  change  in  eccentricity  for  each  orbit 
be  equal  to  zero.  Since  the  initial  orbit  is  circular, 
this  constraint  forces  e  to  remain  zero  for  the  entire  pro¬ 
file,  Radial  thrust  will  be  used  to  negate  changes  in  e 
that  would  be  caused  by  pure  tangential  thrust. 

Although  this  profile  may  not  be  the  absolute  optimum 
because  of  the  additional  constraint,  it  will  be  shown  that 
it  is  at  least  near  optimum. 

The  requirement  that  e  =  0  causes  the  planetary  equa¬ 


tions  to  become 


Fig.  1 .  Shadow  Boundaries 


shadow  exit  point.  Also,  g  will  be  a  measure  of  the  angle 
from  shadow  exit  to  the  next  shadow  entry  and  m  will  be 
one-half  of  the  angle  of  shadow.  See  Figure  1  for  a  deplc 
tion  of  these  angles.  Substituting  s  for  CO  gives 

dfl  W  sin(f+s) 

dt  n  a  sin  i 

di  W  cos(f+s) 

dt  n  a 

Since  n  =  df/dt  ,  a  change  of  independent  variable  can  be 


made  from  t  to  f*  Substitution  yields 


d 0  W  sin(f+s) 

df  n2  a  sin  i 

di  W  cos(f+s) 

df  n2  a 

de  U  sin  f  +  2  V  cos  f 

df  n2  a 

da  2  V 

df  n2 


Substituting  Eq  (6)  gives 


dQ  W  a2  sin(f+s) 


To  determine  the  changes  in  the  orbital  parameters 
for  one  orbit,  these  equations  should  be  integrated  from  0 
27T  •  But  U,  V,  and  W  are  zero  when  the  spacecraft  is 
in  shadow;  so  these  equations  can  be  integrated  from  0  to 
g  where  g  is  the  true  anomaly  at  the  point  where  the  shadow 
is  entered.  The  changes  in  orbital  elements  are  then 


^  fg  W  a2  sin(f*s) 

J 0  M  sin  i 


df 


Ai 


rs  w 

Jo 


g  *  a2  cos(f+s)  df 


i 


g  a2 

Ae  -  /  - (u  sin  f  +  2  V  cos  f)  df 

M 


Aa  = 


i 


g  2  V  a3  df 


The  thrust  accelerations  U,  V  and  W  can  be  modeled  as 
functions  of  f.  From  Figure  2, 


( 


U  =  T  cosy  cos  Oi 
V  =  T  cosy  sin  Q< 


( 


Fig  2.  Acceleration  Components 


where 


7  =  7  it) 
a  =  ou  f) 


Substitution  yields 


T  a2  rs 

Afi - —  sin  7  sin(f+s)  df 

M sin  i  Jo 


cos(f+s)  df 


T  a2  fg 

Ae  - - /  (cosy  cosOi  sin  f 

M  Jn 


*  2  cos^  sinO^  cos  f)  df  (34) 


2  T  a3  rs 

A  a  = - /  cosy  sin  O'  df 

M  In 


(35) 


To  find  the  optimum  thrust  history  for  one  orbit,  the 
approach  used  is  to  maximize  Ai  for  a  given  Aa,  subject  to 
the  additional  constraint  that  Ae  -  0,  The  performance  index 
with  constraint  relationships  is 


j(a 


r 6  T 

,y)  -/  - 

Jo 


g  T  a2 


siny  cos(f*s)  df 


fg  2  T  a3 

Aj  /  - cosy  sinfl^  df  -  Aa 

Jo  M 


Jo 


8ip  a2. 


-(cosy  cosC?  sin  f  +  2cosy  sin 01  cos  f)  df  (36) 


where  aad  X2  are  multipliers.  Simplifying, 
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r  rw  f  * 


J (Oi  ,y  )  = J  J-— ~ — jsin7  cos(f+s)  ♦\1  a  2  cos7  sin# 

♦  X^Ccob'Y  cob  Oi  sin  f  +  2  cos'y  sin#  cos  f)  df 

-  XjAa 


Call  the  integrand  F  for  convenience.  The  calculus  of 
variations  can  be  used  to  show  that  the  above  functional 
has  a  stationary  value  when  the  following  Euler  equations 

Q 

are  satisfied* 


3F  d  [Bf‘ 

a5“ If  ac?’ 


BF  d  IBF 


a7  df lar 


where  primes  indicate  the  derivative  with  respect  to  the  in 
dependent  variable,  f.  Since  Oi'  and  7 '  d0  not  occur  in  F, 
The  Euler  equations  become  algebraic  equations  rather  than 
differential  equations.  Therefore,  J  has  an  extremal  when 


After  performing  the  differentiation, 


3F  \i  2  T  * 


aa  M 


cosy  cos  01 


+ — i - (2  cosTy  cos  Of  cos  f  -  cos^  sinQtf  sin  f) 

M 


3F  T  a2  Xi2  T  a 


-  cosT  cos(f+s)-  — - siny  sinCtf 

ay  M  M 


XaT 


•  (siny  cos  QL  sin  f  +  2  siny  sin  Of  cos  f) 


Substituting  Eqs  (42)  and  (43)  into  Eqs  (40)  and  (41 )  and 
then  simplifying  gives 


tan  d  = 


2(  X^a  +  X2cos 


XpSin  f 


tan  y 


cos(f+s) 


[X|sin2f  ♦  4<  Xja  +  X2C0S  f)2]^ 


Therefore,  the  optimal  control  law  becomes 


Ql  =  tan 


[  X2sin  f  J 

r  cos<f+s) 

1 

a  +  X2coe 


Il'nK'JiM*] 


7)  into  Eqs  (; 


,+  s)  df 


[X|ain2f+4(  Xia+  X2C0S  f)2+cos2(f+s) 
’  cos2(f+s)  df 


[X|ain2f+4(  Xia+  \zcos  f)2+cos2 
[X2sin2f+4-cos  f  (Xia+X2COB  f 


[X|sin2f+4(  Xia+ *  2COS  f)  -*-cos 


( Xia+  \?coe  df 


Fig  3.  Sun-Orbit  Geometrical  Relationship 


As  mentioned  earlier,  s  is  a  phase  angle  measured  in 
the  orbit  plane  between  the  ascending  node  and  the  point 
where  the  spacecraft  exits  the  earth’s  shadow.  Link^  has 
derived  the  following  expressions  for  the  orbital  coordinates 
of  the  anti-sun  point  to  the  orbital  plane. 


sin  B  =  -sin<50  cos  i 
sin  L  =  -cos (0to-  Q) 


♦  cos(50  sin  i  sin(Q^-  Q  ) 

c°s50 

cos  B 


(5- 

(5; 


where  B  and  L  are  angles  as  shown  in  Fig  3,  Q^and  <50  are 
the  right  ascension  and  declination  of  the  sun.  The  same 
reference  also  gives  an  expression  for  half  the  shadow 


Fig  4.  Shadow  Geometry 


angle,  m: 

CO 6  Oo 

cos  m  =  - 

cos  B 

where  Oo  is  half  of  the  maximum  shadow  possible  for  a 
given  altitude. 

Oq  -  7T0*  7Tg  -  R0 

where  7T©  and  7Tg  are  the  parallaxes  of  the  sun  and  sat¬ 
ellite  respectively  and  R  is  the  angular  semi-diameter 
of  the  sun.  See  Figure  4*  It  can  be  seen  in  Figure  3  that 


When  making  the  proceeding  calculation,  care  must  be  exercised 
in  choosing  the  proper  quadrant  for  L.  Once  s  is  known,  the 
upper  limit  of  integration,  g,  can  be  found: 

g  =  27 T-  2  m  =  2( 7T-  m)  (57) 

But  first,  for  convenience,  define 

u  =  a  (58) 

and 

D  =  \|sin2f  ♦  4(u  +  X2cos  f)2+  cos2(f+s)  (59) 

Now  the  optimum  control  history  for  a  single  orbit  can 
be  found  in  the  following  manner.  Given  the  orbital  elements 
and  the  sun-earth  geometrical  relationship,  s  and  g  can  be 
computed  from  Eqs  (56)  and  (57).  Then  the  two  constraint 
relationships, 


g  XpSin2f  +  4  cos  f  (u  +X?cos  f) 


and 


can  be  solved  numerically  to  determine  u  and  X2»  These 
two  values  can  then  be  substituted  into  the  control  law 
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which  i 8  repeated  here; 


Oi  =  tan 


-1 


2  (u  ♦A  2cos  f) 


Apein  f 


(46) 


y  =  tan 


-1 


cos(f+s) 


[A|sin2f  +  4(u  ♦  A2cos  f)2]^ 


(47) 


to  find  the  control  profile.  This  profile  will  maximize 
Ai  for  a  given  Aa,  subject  to  Ae  =  0. 

Examples 

Several  examples  are  presented  here  to  show  the  effects 
caused  by  various  amounts  of  shadow  and  varying  values  of  u. 
In  all  the  examples,  sun-earth  relative  geometry  is  assumed 
such  that 


37T 

s  - - +  m  (62) 

2 


This  value  of  s  corresponds  to  the  following  conditions: 

Qfm*  2?°°  o 
4  =  -23.47° 


a  =  1.03  DU 

n=  iso* 


These  conditions  would  exist  if  the  spacecraft  were  established 
in  a  200  km  parking  orbit  on  the  first  day  of  winter. 


Figure  5  shows  the  effect  of  shadow  width  on  A2» 
Notice  that  when  there  is  no  shadow  (m  =  0),  then  A2  =  0. 
When  A2  s  0  >  the  constraint  equation,  Ae  =  0  ,  is  not 


necessary 


Fig  5.  Shadow  Effects  on  \2 
Figure  6  shows  the  effect  of  shadow  width  on  Aa  and 


Ai,  As  would  be  expected,  decreasing  shadow  width  is 
accompanied  by  increasing  Aa  and  Ai  as  a  result  of  longer 
thrust  application.  Notice  that  u  =  0  corresponds  to 
inclination  change  only  and  that  u  =  oo  corresponds  to 
semi-major  axis  change  only.  For  computational  purposes, 
u  =  100  can  be  used  for  u  =  oo  ,  with  good  numerical 
accuracy,  A  plot  of  A Q  versus  u  is  not  shown  since  aQ 
was  very  near  zero  for  all  cases  that  were  examined.  Also, 
since  aQ  was  always  at  least  three  orders  of  magnitude 
less  than  A  a  and  Ai,  its  equation  will  be  neglected  in 
the  remaining  derivation. 
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A  comparison  of  the  no  shadow  case  and  Alfano's  results 
will  help  verify  the  derivation  in  this  paper.  As  mentioned 
previously,  when  m  =  0  ,  then  \ 2  =  0.  For  \2  =  0  , 

Oi  =7T/2.  and  7  =  tan*1  £cos(f+s)/2uj .  These  values  corre¬ 
spond  exactly  with  those  presented  in  Alfano’s  paper. 
Additionally,  consider  using  the  control  to  cause  semi¬ 
major  axis  change  only.  This  case  corresponds  to  u  =  oo  • 

By  taking  the  limit  as  u— *oo  ,  it  can  be  shown  that 
Aayu./T  a^  =  27T.  Finally,  when  only  change  in  inclination 

p 

is  desired,  Ai/u./  T  a  =  4«  Both  of  the  last  two  results 
also  agree  exactly  with  those  presented  in  Alfano's  paper. 
Therefore,  the  control  law  is  validated  for  the  special 
case  of  no  shadow. 

Figures  7  through  9  show  the  effect  of  the  control,  u, 
for  different  values  of  m,  the  half-shadow  angle.  The  y 
curves  appear  as  would  be  expected;  for  u  =  0,  all  thrust 
is  directed  normal  for  the  orbit  plane;  for  greater  values 
of  u,  less  thrust  is  directed  in  the  normal  direction.  The 
a  curves  may  appear  contradictory  at  first.  As  u  increases 
one  would  expect  more  thrust  to  be  directed  tangentially  to 
achieve  a  greater  Aa.  But  the  figures  indicate  that 
greater  values  of  u  demand  even  greater  divergence  from 
pure  tangential  thrust  ( Oi  =7 T/2).  Recall,  however,  that 
Oi  only  determines  the  direction  of  the  thrust  component 
in  the  orbit  plane. 

The  magnitude  of  that  component  is  T  cos')'  •  As  u 
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increases,  cosy  increases*  So  even  though  it  appears  that 
the  tangential  component  decreases  due  to  Q  ,  the  coupling 
through  the  7  dependence  insures  that  the  tangential  force 
is  increasing*  Of  course,  the  increasing  tangential  accel¬ 
eration  requires  more  radial  accelerations  to  maintain 
Ae  =  0;  thus  causing  it  to  appear  that  less  tangential  force 
is  being  applied. 


Ill  Slow  Timescale  Problem 


Problem  Statement 

The  purpose  of  solving  the  slow  timescale  problem  is  to 
determine  how  much  to  change  semi -major  axis  and  inclination 
on  each  orbit  so  that  final  boundary  conditions  are  reached 
in  minimum  time*  Fast  timescale  results  are  used  to  ensure 
optimal  control  during  each  orbit  and  to  define  the  amount 
of  change  possible  on  a  given  orbit.  Mass  is  recomputed 
for  each  orbit  to  compensate  for  propellant  loss. 

Derivation 

Before  the  minimum  time  control  problem  can  be  solved, 
a  few  preliminary  steps  will  be  taken.  An  expression  is 
needed  for  da/dt  and  di/dt  for  the  slow  timescale  problem. 

Since  the  orbital  elements  change  so  little  on  each  orbit, 
these  rates  can  be  approximated  by 

da  Aa 

- ^ -  (63) 

dt  At 

di  Ai 

- » -  (64) 

dt  At 


where  At  is  the  elapsed  time  for  the  particular  orbit.  For 
a  circular  orbit, 


At 


TP* 


(65) 


But  the  presence  of  shadow  makes  the  above  approximations 


inaccurate.  The  elements  a  and  i  will  only  change  when  the 
spacecraft  is  in  sunlight.  Therefore,  At  must  be  adjusted 
to  only  cover  that  portion  of  the  orbit  when  thrust  is 


being  applied. 


At  - - TP=  2(7 T-  m) 

27 T 


V| 
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Substitution  yields 


da  2  T  a^2  fs  (u  +  X?cos  f )  df 


dt  /u2(7T-m) 


cos  (f+s)  df 


dt  2/lx  (7T-m)J0 


Since  the  problem  is  to  find  a  minimum  thrusting  time  and 
hence  minimum  fuel,  successive  orbits  will  combine  as  if 
there  were  no  coasting  time  through  the  shadow.  The  coast 
ing  time  will  be  computed  so  the  total  transfer  elapsed 
time  will  be  known  at  the  end. 

As  mentioned  previously,  corrections  will  be  made  for 
each  orbit  to  compensate  for  increased  acceleration  caused 
by  decreasing  mass.  For  a  constant  thrust  ion  rocket,  pro 
pellant  flow  is  constant.  Acceleration  as  a  function  of 
time  can  be  modeled  as 


hnl 


where  T0  is  the  initial  vehicle  acceleration,  t  is  the  time 
and  Ap  is  the  specific  mass  flow  rate  (mass  flow  rate  divided 
by  initial  vehicle  mass). 

Since  Eqs  (67)  and  (68)  can  be  solved  more  easily  if 
there  is  no  time  dependence,  a  change  of  variable  can  be 
made.  A  new  independent  variable  T  can  be  defined  such  that 


dT 


Integrating,  with  T  =  0  when  t  =  0  yields 


(7< 


T 

T  =  -  In  (1  -  m  t)  (7 

ap 

T is  the  total  accumulated  velocity  change.  Minimizing 
T  will  minimize  thrusting  time  and  fuel  expended. 

Converting  from  t  to  T  gives 


g  (u  +  \2cos  f)  df 

di  a^  re  cos2(f+s) 
dT  /jl  ^2(7T-m) JQ 

Bryson  and  Ho10  have  shown  that  a  minimum  time 


da  2  a3/2  C 
dT  /x*(7T-m) J0 


(7« 


(7; 


solution  satisfies  the  following  conditions: 


=  0 


^  =  1»2,3»»»*,r 


(75) 


3H 

=  al 


j  =  1,2,3»***,q 


(76) 


where  H  is  the  Hamiltonian,  X-;'8  are  Lagrange  multipliers, 
Xj*8  are  the  state  variables,  and  the  u^’s  are  the  control 
variables*  In  this  problem,  r  =  1  producing  only  one 
optimality  condition  and  q  =2  because  the  simplified 
system  has  only  two  degrees  of  freedom.  These  equations 
will  now  be  applied  to  this  specific  problem. 

For  this  problem,  the  Hamiltonian  is 


da  di 

H  *  ' 


(77) 


Since  t  does  not  appear  explicitly  in  H, 


H  =  0 


(78) 


Therefore,  H(tf)  =  0  implies 


H(t)  =  0 


(79) 


for  all  t  2  0  .  So 


da  di 

Aa—  +  X±— r  +  1  =  0 
a  d T  1  dT 


for  all  time.  Substitution  produces 
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s  [X?sin2f  +  cos2(f+s)l  df 


a i  iii 


UT(7T-m) 


Xpcos  f)cos2(f+s)  df 


For  this  problem,  after  the  change  of  independent  variable, 
Eq  (76)  becomes 

.  .  3H 

Ao  -  ^ 


where  the  prime  indicates  differentiation  with  respect 
to  T.  Now 


3li  ^  a_N,  ^ 

3a  "  Aa  9a  ST  1  3a  dT 


But,  remembering  that  s  is  a  function  of  m  and  m  is  a  function 
of  both  a  and  i, 


3  da 
Ba  dT 


f~B  da]  3  dal  3  m 

Ba  d7l  Bm  dT  Ba 


(90) 


3  da  3  f  s  (u  +  X2C0S  f )  df 

3a  dT  JQ  D^~ 


Since  g  is  also  a  function  of  m 


2  a^2(u  +  Xpcos  f)  df 


3  da  rg  3  2a3/2(u*Xi 

3m  dT  Jq  3m  fJL^(TT-m) 


2  a5//2(u  +  X2C0S  f) 

yU  ^  ( 7T  -m) 


6 


3S 

3m 


(91) 


And,  finally 


3  da  2  a 


3/2 


9a  dT  "  ^(.TT-mV 


r 

Jo 


(u  +  X?cos  df 


2  a 


^/^  r&  (u  +  Xpcos  f )sin(f+s)cos(f+s)  df 


— r 

(7T-m)  JQ 


4  a^2(u+  Xpcos  2m) 


7T-m)[X2sin22m+4(u+  X2cos2m)2+cos2(s-2m)] 


(92) 


Now,  to  find  3m/ 3a,  recall  Eqs  (52),  (54) »  and  (55)* 
Substitution  of  the  appropriate  physical  constraints  gives 

cosO^  = 


1  -■ 


Fo  +  “  G0 
a 


(93) 
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where  F0  and  G0are  constants  depending  on  the  size  of  the 
earth  and  sun  and  the  distance  between  them. 

F0  =  .99994668  G0  =  .0046229536 

Also, 

cos  B  =  £l  -(cos(5#  sin  i  sin(O^-Q)-  sin<5^  cos  i)2]^  (94) 


B  di  1 

3a  d T  4  /uL^a^i  7T -m) 


.g  2 

o  r»n  o 


cosc(f+s)  df 


(99) 


3  di_  -a* 

3m  d T  "  /ut*(7T-m) 


sin(f+s)  cos(f+s)  df 


a^  cos2(f+s)  df 

+  2M*(7T-»)2  JQ  D* 


a^  fs  sin(f+s)  cos^(f+s)  df 


2M*(7T- m) 


=rjc 


■^72 


a*cos  (s-2m) 


fj^(7T-m)  [X|sin22m+4(u+  X2cos2m)2+cos2(s-2m)]^ 


(100) 
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Now,  Eqs  (81 )  and  (85)  can  be  solved  simultaneously 
for  X  5111(1  X  •  in  terms  of  integrals  containing  the  state 

SI  X 

and  control  variables.  From  Eq  (85), 


Substituting  this  result  into  Eq  (81)  gives 


These  two  equations  can  be  substituted  into  either  equation 
(101)  or  (102)  to  also  produce  an  integral  equation  in  terms 
of  the  state  and  control  variables.  Now,  using  the  new 
expression  for  Eqs  (101)  or  (102),  and  the  two  state  Eqs 
(72)  and  (73),  the  minimum  time  problem  can  be  solved. 


There  are  three  equations  and  three  unknowns:  a,  i,  and  u. 
Solving  these  equations  is  not  a  trivial  task,  however. 
Implementation  of  this  solution  is  described  in  the  following 
section. 


Implementation 

The  minimum  time  profile  is  determined  in  the  following 
manner.  Choose  Eq  (101)  or  (102)  to  use.  To  demonstrate, 
use  Eq  (101).  Pick  a  starting  value  for  X  •  Given  this 

9i 

X  _(0),  Eq  (105)  can  be  solved  numerically  to  determine  what 

a 

u  would  produce  that  X»«  Remember  that  selecting  or  find- 
ing  u  also  determines  X2  as  a  result  of  satisfying  the  con¬ 
straint  Ae  =  0  .  Now  these  values  of  u  and  X2  are  used  in 
the  state  Eqs  (72)  and  (73)  to  find  Aa  and  Ai.  Also,  a,  i 
and  m  are  used  to  compute  At  and  AT.  This  AT  is  used  to 
find  the  new  Xa«  Xa  is  then  used  to  repeat  the  process 
until  final  boundary  conditions  are  met.  If  final  boundary 
conditions  are  not  met,  a  new  Xo(0)  is  chosen.  This  pro- 

SL 

cess  continues  until  a  X.(°)  is  found  which  causes  the 
final  boundary  conditions  to  be  achieved.  Then  u  will  be 
determined  for  the  entire  profile. 

Although  this  method  of  solution  may  sound  rather 
straight  forward,  there  can  be  some  major  problems.  At 
first  glance,  it  would  appear  that  one  should  choose  the  Xi 
equation  to  use.  That  equation  has  fewer  terms  and  Y  as  a 
factor  in  all  its  terms.  Recall  that  Y  is  the  partial  of  ra 
with  respect  to  i.  When  ra  =  0,  Y  =  0  and 


(107) 


Therefore 

X  ^  =  constant  (108) 

It  was  found  by  experience,  however,  that  in  the  range  of 
interest  for  values  of  \^,  u  is  double-valued.  It  was 
also  found  that  \ ^  varied  dramatically  for  different  values 
of  a  and  i.  See  figures  10  and  11.  Figure  10  shows  X_ 
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and  \±  versus  u  for  a  =  1.03  Oil’s  and  i  =  .4974  radians. 

Figure  11  shows  the  same  functions  for  a  =  6.6  DU's  (geo¬ 
synchronous)  and  i  =  0.  Notice  that  the  X  curve  main- 
tains  its  shape  and  only  changes  slightly.  The  X^  curve, 
however,  is  not  as  well  behaved.  The  peak  on  that  curve 
shifts  upward  and  to  the  right.  Using  X^»  the  optimal 
profile  attempts  to  move  from  the  right  side  of  the  peak 
to  the  left.  Numerical  search  schemes  in  the  vicinity  of 
the  zero  slope  diverged  and  therefore  an  optimal  profile 
could  not  be  found. 

By  using  the  Xa  equation,  an  optimal  profile  could 
be  found  as  long  as  large  values  of  u  are  not  encountered. 

This  would  occur  if  the  starting  Xa  was  chosen  where  the 
Xa  curve  was  nearly  horizontal.  Fortunately,  that  region 
corresponds  to  transfers  which  will  be  shown  to  not  be 
optimal. 

The  second  major  problem  in  implementing  this  solu¬ 
tion  is  developing  a  search  scheme  to  find  the  correct 
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Fig  10.  \a  8111(1  Xi  *or  Initial  Conditions 
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\  (0).  For  all  cases  considered  in  this  paper,  choosing 
a  ^a(0)  near  zero  caused  i^  to  be  reached  before  a^,.  For 
\a(0)  in  the  horizontal  portion  of  its  graph,  af  was 
reached  first*  The  optimum  choice  of  for  a  minimum 

time  solution  occurs  somewhere  between  these  other  values. 
Unfortunately,  more  than  one  local  minimum  may  exist  in 
this  region.  In  the  next  chapter,  it  will  be  shown  that 
this  does  not  appear  to  be  a  serious  problem  for  the  par¬ 
ticular  transfers  that  were  considered.  When  more  than 
one  minimum  occured,  the  transfer  times  were  very  close 
to  one  another. 

In  summary,  the  method  outlined  in  the  chapter  can 
be  used  to  solve  the  minimum  time  transfer  problem.  The 
solution  may  only  be  a  local  minimum,  however. 


IV  Results 


The  equations  derived  in  the  proceeding  chapters  were 
applied  to  three  specific  transfers.  One  transfer  did  not 
involve  any  time  in  shadow;  another  caused  the  spacecraft 
to  be  in  shadow  for  about  half  its  orbits;  and  the  last 
transfer  required  shadow  penetration  during  about  95%  of 
its  orbits.  All  three  transfers  required  the  same  semi- 
major  axis  and  inclination  change.  For  all,  ai  correspond¬ 
ed  to  a  parking  orbit  altitude  of  200  km  and  af  correspond¬ 
ed  to  geo-synchronous  altitude.  Beginning  and  ending  in¬ 
clinations  were  chosen  to  cause  the  varying  amounts  of 
shadow  time.  In  all  three  cases, 

do-  270*  n  =  ,80’  <5o=  -23.47" 

Also,  to  aid  comparison,  Alfano's  values  for  specific  im¬ 
pulse  (5000  sec)  and  specific  mass  flow  rate  (2.0  x  10-^/8ec) 
were  used.  Figures  12  through  14  show  how  a,  i  and  shadow 
angle  vary  during  each  of  the  three  profiles. 

Case  I 

This  case  was  the  transfer  that  did  not  penetrate  the 
earth' 8  shadow  at  any  time.  This  case  was  designed  primar¬ 
ily  to  validate  the  slow  timescale  solution  by  comparing 
it  with  Alfano's  results  for  the  identical  transfer.  Using 
the  solution  presented  here,  the  total  accumulated  velocity 
change  (T)  was  only  2.7%  greater  than  Alfano's.  The  differ¬ 
ence  was  caused  primarily  by  numerical  problems  associated 
with  finding  A 2  from  the  constraint  relationship,  Ae  =  0. 
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Fig  15.  Case  II  Initial  Conditions 


Although  \2  should  equal  zero  when  there  is  no  shadow,  X2 
was  always  found  to  be  very  near  zero.  Consequently,  a  small 
amount  of  radial  thrust  was  applied  throughout  the  transfer 
making  it  less  efficient.  A  more  accurate  search  routine 
for  X2  should  reduce  or  eliminate  the  difference.  For  this 
reason,  it  is  believed  that  this  solution  is  validated. 

Case  II 

The  geometry  for  this  case  is  shown  in  Figure  15.  The 
initial  conditions  would  occur  with  a  noon  launch  out  of 
Cape  Canaveral  on  the  first  day  of  winter.  In  this  transfer, 
the  spacecraft  experiences  an  eclipse  on  each  of  its  first 
20  orbits. 


Fig  16.  Case  III  Initial  Conditions 

For  this  transfer,  T  was  only  1.17%  greater  than  the 
case  with  no  shadow.  This  small  amount  represents  the 
penalty  associated  with  maintaining  ^e  =  0  during  the  orbits 
when  the  shadow  is  penetrated.  Even  the  total  transfer  time 
including  coasting  periods  is  only  1*0  TU's  or  5.62%  longer. 

Case  III 

The  initial  geometry  for  this  transfer  is  shown  on  Fig¬ 
ure  16.  These  conditions  would  occur  with  a  midnight  launch 
on  the  first  day  of  winter.  This  transfer  was  chosen  be¬ 
cause  it  required  the  spacecraft  to  penetrate  the  shadow 
for  many  orbits.  As  expected,  this  transfer  required  more 
fuel  and  time.  There  is  a  3.76%  increase  in  T.  The  total 
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time  experienced  a  large  24*8%  increase,  however.  Again, 
this  is  not  surprising,  since  the  spacecraft  had  many  more 
coasting  periods. 

Multiple  Local  Minima 

As  previously  discussed,  several  minima  may  exist. 

Using  Case  II  as  an  example,  it  was  found  that  at  least  two 
local  minima  did  exist.  For  \_(0)  =  -.69023,  T  and  thrusting 

A 

time  were  as  described  previously  and  shown  in  Figure  12. 

For  \a(0)  =  -*73123,  T  was  only  ,001  greater.  The  two 
profiles  were  slightly  different.  In  the  second  profile, 
semi-major  axis  changed  more  rapidly  at  first  causing 
shadow  exit  one  orbit  earlier. 

The  presence  of  multiple  minima  is  not  really  a  pro¬ 
blem,  however.  A  reasonable  step  size  for  \  (0)  of  .01 
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identified  where  the  minima  could  be  found.  Also,  it  was 
noted  that  there  is  a  specific  range  of  \  (0)  where 
these  minima  seem  to  occur.  If  Xa^0^  is  such  that  if 
is  reached  first,  and,  if  thrust  is  used  after  that  point 
to  change  a  only,  then  the  amount  of  time  needed  to  reach 
is  an  indication  of  whether  a  minimum  can  exist  in  a 
neighborhood  of  that  Xft(0)*  If  af  cannot  be  reached 
within  one  orbit,  then  \a(0 )  does  not  appear  to  be  near 
a  minimum.  A  plot  of  T  versus  \a(0)  seems  to  be  steadily 

A 

decreasing  in  this  region;  see  Figure  17.  The  same  argu¬ 
ment  can  be  made  for  XQ(°)'S  where  a„  is  reached  first. 

If  if  cannot  be  reached  within  one  orbit  using  only  thrust 


Fig  17.  T  versus  \  (0) 
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normal  to  the  plane,  then  a  minimum  does  not  apparently 
exist  in  that  region.  Again,  T seems  to  be  steadily  in¬ 
creasing  in  this  area.  In  either  case,  if  the  other 
boundary  condition  can  be  reached  within  one  orbit,  then 
a  minimum  is  nearby.  Also,  the  new  \_(0)  should  be 
greater  if  if  was  reached  first  and  less  if  was 
achieved  first.  Knowing  these  characteristics  makes  it 
relatively  easy  to  find  the  minima. 

General  Observations 

There  are  two  additional  characteristics  of  these 
transfers  that  bear  mentioning.  Several  other  initial 
conditions  were  used  to  examine  certain  characteristics. 


Full  transfers  were  not  accomplished,  but  two  trends  did 
appear. 

The  first  trend  to  note  is  that  in  no  case  did  the 
control  law  direct  an  inclination  change  away  from  i^. 

It  was  thought  that  the  minimum  time  solution  may  require 
a  Ai  away  from  i^  to  reduce  shadow  time.  Apparently,  the 
cost  of  moving  away  from  the  final  boundary  condition  is 
more  than  the  benefit  gained  from  reducing  the  shadow. 

The  other  trend  that  was  noticed  was  that  the  control, 
u,  increases  until  the  shadow  is  less  than  approximately 
2.1  radians.  If  the  shadow  is  less  than  that  amount,  u 
decreases.  The  value  at  which  this  occurs  is  not  constant 
for  all  the  transfers,  but  it  is  near  2.1  radians. 

It  should  also  be  pointed  out  that  these  trends  were  not 
contradicted  in  any  of  these  trials. 
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Conclusions  and  Recommendations 

A  minimum  time,  minimum  fuel  control  law  for  a  vehicle 
using  discontinuous  low  thrust  has  been  presented.  This 
control  is  optimal,  subject  to  the  constraint  that  a®  =  0. 

•-If  vehicle  design  prevents  continuous  thrust  applica¬ 
tion,  there  is  only  a  small  penalty  over  a  design  which 
allows  continuous  thrust.  In  particular,  an  ion  rocket 
powered  by  solar  panels  would  be  competitive  with  one 
powered  by  nuclear  power,  even  if  the  former  could  not 
operate  in  shadow.  The  reduction  in  weight  and  complexity 
by  eliminating  a  nuclear  power  source  may  more  than  offset 
the  weight  of  a  small  amount  of  additional  fuel. 

This  solution  assumed  a  spherical  earth,  as  well  as 
constant  Q^and  <50.  The  oblateness  of  the  earth  would 
cause  the  line  of  nodes  to  regress.  Additionally,  the 
actual  changes  in  0io  and  do  will  affect  the  size  of  the 
shadow.  These  two  effects  combine  to  make  the  length  of 
the  transfer  highly  dependent  on  the  launch  time  and  in¬ 
clination.  The  proper  selection  of  when  to  launch  into 
which  initial  parking  orbit  can  minimize  the  effect  of 
shadow  penetration. 

This  solution  also  assumed  low  thrust.  If  thrust 
is  too  great,  this  control  law  is  invalid  because  Aa  and 
Ai  would  not  be  small  on  each  orbit. 


A  recommendation  to  further  this  study  would  be  to  find 
an  optimum  launch  time  to  perform  a  given  transfer*  Another 
continuation  would  be  to  find  the  optimum  control  law  with¬ 
out  requiring  Ae  =  0. 
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